Roider DLBCL scRNA-seq reprocessing

Rendered 12/12/2025

Author

Jeremy M. Simon

Set up Seurat workspace

# Load libraries
library(data.table)
library(devtools)
library(presto) 
library(glmGamPoi)
library(sctransform)
library(Seurat) 
library(tidyverse)
library(miQC)       
library(SeuratWrappers) 
library(flexmix)
library(SingleCellExperiment)
library(SummarizedExperiment)
library(readxl)
library(Matrix)
library(speckle)
library(scater)
library(patchwork)
library(vctrs)
library(harmony)
library(scDblFinder)
library(ggridges)
library(scGate)
library(ggsankey)
library(ggpubr)

# Set global options for Seurat v5 objects
options(Seurat.object.assay.version = 'v5')

Import 10x data and create Seurat objects with sample names prepended to cell barcodes

DLBCL1 <- Read10X("DLBCL1",strip.suffix = T)
DLBCL2 <- Read10X("DLBCL2",strip.suffix = T)
DLBCL3 <- Read10X("DLBCL3",strip.suffix = T)
FL1 <- Read10X("FL1",strip.suffix = T)
FL2 <- Read10X("FL2",strip.suffix = T)
FL3 <- Read10X("FL3",strip.suffix = T)
FL4 <- Read10X("FL4",strip.suffix = T)
rLN1 <- Read10X("rLN1",strip.suffix = T)
rLN2 <- Read10X("rLN2",strip.suffix = T)
rLN3 <- Read10X("rLN3",strip.suffix = T)
tFL1 <- Read10X("tFL1",strip.suffix = T)
tFL2 <- Read10X("tFL2",strip.suffix = T)

DLBCL1.seurat <- CreateSeuratObject(counts = DLBCL1)
Warning: Feature names cannot have underscores ('_'), replacing with dashes
('-')
DLBCL1.seurat <- RenameCells(DLBCL1.seurat,add.cell.id = "DLBCL1")
DLBCL2.seurat <- CreateSeuratObject(counts = DLBCL2)
Warning: Feature names cannot have underscores ('_'), replacing with dashes
('-')
DLBCL2.seurat <- RenameCells(DLBCL2.seurat,add.cell.id = "DLBCL2")
DLBCL3.seurat <- CreateSeuratObject(counts = DLBCL3)
DLBCL3.seurat <- RenameCells(DLBCL3.seurat,add.cell.id = "DLBCL3")
FL1.seurat <- CreateSeuratObject(counts = FL1)
Warning: Feature names cannot have underscores ('_'), replacing with dashes
('-')
FL1.seurat <- RenameCells(FL1.seurat,add.cell.id = "FL1")
FL2.seurat <- CreateSeuratObject(counts = FL2)
Warning: Feature names cannot have underscores ('_'), replacing with dashes
('-')
FL2.seurat <- RenameCells(FL2.seurat,add.cell.id = "FL2")
FL3.seurat <- CreateSeuratObject(counts = FL3)
Warning: Feature names cannot have underscores ('_'), replacing with dashes
('-')
FL3.seurat <- RenameCells(FL3.seurat,add.cell.id = "FL3")
FL4.seurat <- CreateSeuratObject(counts = FL4)
FL4.seurat <- RenameCells(FL4.seurat,add.cell.id = "FL4")
rLN1.seurat <- CreateSeuratObject(counts = rLN1)
Warning: Feature names cannot have underscores ('_'), replacing with dashes
('-')
rLN1.seurat <- RenameCells(rLN1.seurat,add.cell.id = "rLN1")
rLN2.seurat <- CreateSeuratObject(counts = rLN2)
rLN2.seurat <- RenameCells(rLN2.seurat,add.cell.id = "rLN2")
rLN3.seurat <- CreateSeuratObject(counts = rLN3)
rLN3.seurat <- RenameCells(rLN3.seurat,add.cell.id = "rLN3")
tFL1.seurat <- CreateSeuratObject(counts = tFL1)
Warning: Feature names cannot have underscores ('_'), replacing with dashes
('-')
tFL1.seurat <- RenameCells(tFL1.seurat,add.cell.id = "tFL1")
tFL2.seurat <- CreateSeuratObject(counts = tFL2)
tFL2.seurat <- RenameCells(tFL2.seurat,add.cell.id = "tFL2")

Retrieve published cell-cluster metadata

This is seemingly on a downsampled set of cells, perhaps done so for easier plotting in the shiny app

download.file("https://github.com/anders-biostat/single-cell-lymphoma-explorer/raw/master/metaInfo.rds", 
              destfile = "Roider_DLBCL_metaInfo.rds")
roider.meta <- readRDS("Roider_DLBCL_metaInfo.rds")

head(roider.meta$bcells)
                          UMAP1     UMAP2 Subset
tFL2_GATTCAGAGCTGGAAC -3.726201 -7.181145      5
tFL2_AACACGTCACACAGAG -3.260301 -5.485855      5
tFL2_ACATCAGGTTCCGGCA -4.213102 -6.762125      5
tFL2_CCCAGTTTCACCGGGT -3.345386 -5.145964      5
tFL2_ACTTGTTTCGTGGACC -4.179855 -6.458187      5
tFL2_CTCGTCATCACCGGGT -2.707596 -6.752877      5

Merge objects

merged.roider <- merge(x = DLBCL1.seurat, y=c(DLBCL2.seurat, DLBCL3.seurat, FL1.seurat, FL2.seurat, FL3.seurat, FL4.seurat, rLN1.seurat, rLN2.seurat, rLN3.seurat, tFL1.seurat, tFL2.seurat))
dim(merged.roider)
[1] 45068 41786

QC filter

merged.roider <- subset(merged.roider, subset = nFeature_RNA > 250)
dim(merged.roider)
[1] 45068 41645
merged.roider <- PercentageFeatureSet(merged.roider, pattern = "^MT-", col.name = "percent.mt")
merged.roider <- RunMiQC(merged.roider, 
                        percent.mt = "percent.mt", 
                        nFeature_RNA = "nFeature_RNA", 
                        posterior.cutoff = 0.7, 
                        model.slot = "flexmix_model")
Warning: Adding a command log without an assay associated with it
merged.roider <- subset(merged.roider, miQC.keep == "keep")
dim(merged.roider)
[1] 45068 40073
data.frame(table(str_replace_all(colnames(merged.roider),"_.+","")))
     Var1 Freq
1  DLBCL1 4008
2  DLBCL2 4927
3  DLBCL3 2412
4     FL1 3668
5     FL2 5090
6     FL3 4686
7     FL4 3147
8    rLN1 3082
9    rLN2 2028
10   rLN3 2779
11   tFL1 1998
12   tFL2 2248

Add meta data

disease <- str_replace_all(colnames(merged.roider),"\\d_.+","")
sample <- str_replace_all(colnames(merged.roider),"_.+","")

merged.roider <- AddMetaData(merged.roider,disease,col.name="Disease")
merged.roider <- AddMetaData(merged.roider,sample,col.name="Sample")

Join then re-split RNA counts layers by Sample

merged.roider[['RNA']] <- JoinLayers(merged.roider[['RNA']])
merged.roider[["RNA"]] <- split(merged.roider[["RNA"]], f = merged.roider$Sample)

Normalize and scale data

Regress out mitochondrial contribution

merged.roider <- NormalizeData(merged.roider)
merged.roider <- FindVariableFeatures(merged.roider, 
                                    assay="RNA", 
                                    layer="data", 
                                    selection.method = "vst", 
                                    nfeatures = 5000)
merged.roider <- ScaleData(merged.roider, vars.to.regress = "percent.mt")

Run PCA

merged.roider <- RunPCA(merged.roider, npcs = 200, verbose = FALSE)
ElbowPlot(merged.roider, ndims = 200, reduction = "pca")

Plot unintegrated UMAP

merged.roider <- RunUMAP(merged.roider, 
                        reduction = "pca", 
                        dims = 1:30, 
                        reduction.name = "umap.unintegrated")
Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
This message will be shown once per session
10:55:58 UMAP embedding parameters a = 0.9922 b = 1.112
Found more than one class "dist" in cache; using the first, from namespace 'BiocGenerics'
Also defined by 'spam'
10:55:58 Read 40073 rows and found 30 numeric columns
10:55:58 Using Annoy for neighbor search, n_neighbors = 30
Found more than one class "dist" in cache; using the first, from namespace 'BiocGenerics'
Also defined by 'spam'
10:55:58 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
10:56:09 Writing NN index file to temp file /tmp/RtmpGN9GPk/file150e28686366a7
10:56:09 Searching Annoy index using 1 thread, search_k = 3000
10:56:25 Annoy recall = 100%
10:56:26 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
10:56:30 Initializing from normalized Laplacian + noise (using RSpectra)
10:56:41 Commencing optimization for 200 epochs, with 1735844 positive edges
10:57:12 Optimization finished
DimPlot(merged.roider, reduction = "umap.unintegrated", group.by = c("Disease", "Sample"))

Call preliminary clusters for the purposes of doublet discrimination

merged.roider <- FindNeighbors(merged.roider, dims = 1:30, reduction = "pca")
Computing nearest neighbor graph
Computing SNN
merged.roider <- FindClusters(merged.roider, 
                             resolution = 0.2, 
                             algorithm = 2)
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 40073
Number of edges: 1491903

Running Louvain algorithm with multilevel refinement...
Maximum modularity in 10 random starts: 0.9798
Number of communities: 21
Elapsed time: 10 seconds
DimPlot(merged.roider, reduction = "umap.unintegrated", label = T)

Identify and remove doublets

This uses raw counts as input

Convert seurat to sce and check colData

# Join RNA layers
merged.roider[['RNA']] <- JoinLayers(merged.roider[['RNA']])

merged.roider.sce <- as.SingleCellExperiment(merged.roider, assay = "RNA")
merged.roider.sce
class: SingleCellExperiment 
dim: 45068 40073 
metadata(0):
assays(2): counts logcounts
rownames(45068): RP11-34P13.3 FAM138A ... AC145212.1 MAFIP
rowData names(0):
colnames(40073): DLBCL1_AAACCTGAGCCATCGC DLBCL1_AAACCTGAGCCTTGAT ...
  tFL2_TTTGTCATCAGGCGAA tFL2_TTTGTCATCAGGTAAA
colData names(11): orig.ident nCount_RNA ... seurat_clusters ident
reducedDimNames(2): PCA UMAP.UNINTEGRATED
mainExpName: RNA
altExpNames(0):
colData(merged.roider.sce)
DataFrame with 40073 rows and 11 columns
                           orig.ident nCount_RNA nFeature_RNA percent.mt
                          <character>  <numeric>    <integer>  <numeric>
DLBCL1_AAACCTGAGCCATCGC SeuratProject       9180         2007    4.65142
DLBCL1_AAACCTGAGCCTTGAT SeuratProject      16902         3280    1.36670
DLBCL1_AAACCTGAGTTCGCAT SeuratProject       4800         1490    5.06250
DLBCL1_AAACCTGCACCACCAG SeuratProject      17052         4003    2.18743
DLBCL1_AAACCTGGTTAAGATG SeuratProject       8116         2297    1.88517
...                               ...        ...          ...        ...
tFL2_TTTGTCAGTCATCGGC   SeuratProject       7075         2077   2.388693
tFL2_TTTGTCAGTTCACGGC   SeuratProject       5901         1677   1.745467
tFL2_TTTGTCATCAACGCTA   SeuratProject       6722         2240   2.424874
tFL2_TTTGTCATCAGGCGAA   SeuratProject       4281         1714   0.981079
tFL2_TTTGTCATCAGGTAAA   SeuratProject       5842         1760   1.763095
                        miQC.probability   miQC.keep     Disease      Sample
                               <numeric> <character> <character> <character>
DLBCL1_AAACCTGAGCCATCGC       0.05811154        keep       DLBCL      DLBCL1
DLBCL1_AAACCTGAGCCTTGAT       0.01431377        keep       DLBCL      DLBCL1
DLBCL1_AAACCTGAGTTCGCAT       0.16824566        keep       DLBCL      DLBCL1
DLBCL1_AAACCTGCACCACCAG       0.00701990        keep       DLBCL      DLBCL1
DLBCL1_AAACCTGGTTAAGATG       0.00628381        keep       DLBCL      DLBCL1
...                                  ...         ...         ...         ...
tFL2_TTTGTCAGTCATCGGC         0.00509071        keep         tFL        tFL2
tFL2_TTTGTCAGTTCACGGC         0.00515423        keep         tFL        tFL2
tFL2_TTTGTCATCAACGCTA         0.00530276        keep         tFL        tFL2
tFL2_TTTGTCATCAGGCGAA         0.01042485        keep         tFL        tFL2
tFL2_TTTGTCATCAGGTAAA         0.00532506        keep         tFL        tFL2
                        RNA_snn_res.0.2 seurat_clusters    ident
                               <factor>        <factor> <factor>
DLBCL1_AAACCTGAGCCATCGC               1               1        1
DLBCL1_AAACCTGAGCCTTGAT               1               1        1
DLBCL1_AAACCTGAGTTCGCAT               1               1        1
DLBCL1_AAACCTGCACCACCAG               1               1        1
DLBCL1_AAACCTGGTTAAGATG               1               1        1
...                                 ...             ...      ...
tFL2_TTTGTCAGTCATCGGC                10              10       10
tFL2_TTTGTCAGTTCACGGC                10              10       10
tFL2_TTTGTCATCAACGCTA                10              10       10
tFL2_TTTGTCATCAGGCGAA                5               5        5 
tFL2_TTTGTCATCAGGTAAA                10              10       10

Run scDblFinder

Set the dbr.sd very high to better allow thresholds to be set based on misclassification rates per sample

merged.roider.sce <- scDblFinder(merged.roider.sce,
                                 samples = "Sample",
                                 dbr.sd = 1,
                                 clusters = "seurat_clusters")

Inspect results

# Look at the classes
table(merged.roider.sce$seurat_clusters, merged.roider.sce$scDblFinder.class)
    
     singlet doublet
  0     6291     215
  1     3932      62
  2     3945      21
  3     3366      82
  4     3371      36
  5     2866      75
  6     2609      16
  7     2561       9
  8     2243      38
  9     2163      13
  10    1635      33
  11     767      76
  12     689     151
  13     631     205
  14     740      62
  15     272     217
  16      97     110
  17     140      39
  18       0     120
  19      88      25
  20      39      23
table(merged.roider.sce$Sample, merged.roider.sce$scDblFinder.class)
        
         singlet doublet
  DLBCL1    3932      76
  DLBCL2    4662     265
  DLBCL3    2354      58
  FL1       3532     136
  FL2       4951     139
  FL3       4434     252
  FL4       3039     108
  rLN1      2833     249
  rLN2      1957      71
  rLN3      2710      69
  tFL1      1874     124
  tFL2      2167      81
# Look at the scores
summary(merged.roider.sce$scDblFinder.score)
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
0.0000757 0.0020391 0.0033650 0.0572434 0.0138431 0.9999380 

Save doublet classifications into main Seurat object

merged.roider$doublet_classification <- merged.roider.sce$scDblFinder.class

Count singlets and doublets

table(merged.roider$doublet_classification)

singlet doublet 
  38445    1628 
table(merged.roider$doublet_classification, merged.roider$seurat_clusters)
         
             0    1    2    3    4    5    6    7    8    9   10   11   12   13
  singlet 6291 3932 3945 3366 3371 2866 2609 2561 2243 2163 1635  767  689  631
  doublet  215   62   21   82   36   75   16    9   38   13   33   76  151  205
         
            14   15   16   17   18   19   20
  singlet  740  272   97  140    0   88   39
  doublet   62  217  110   39  120   25   23

Plot singlets/doublets in UMAP space

DimPlot(merged.roider,reduction = "umap.unintegrated", group.by = "doublet_classification")

Subset object to remove doublets and count remaining cells

merged.roider.singlets <- subset(merged.roider, doublet_classification %in% c("singlet"))
dim(merged.roider.singlets)
[1] 45068 38445
# Count remaining cells per initial cluster
table(merged.roider.singlets$seurat_clusters)

   0    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15 
6291 3932 3945 3366 3371 2866 2609 2561 2243 2163 1635  767  689  631  740  272 
  16   17   18   19   20 
  97  140    0   88   39 

Set up gating strategy to focus only on B-cell lineage

We gate to remove: * CD3+ (T cells) * DERL3+/ITM2C+ (Plasma cells)

# Plot genes to be used in gating model
FeaturePlot(merged.roider.singlets, 
            features=c("CD79A", "CD79B", 
                      "MS4A1", "BANK1", 
                      "CD3E", "CD3D", 
                      "CD2", "GZMK", 
                      "DERL3", "ITM2C"), 
            order = T, 
            ncol = 5)

Set up gating model and run scGate

scGate_model <- gating_model(name = "Tcells", 
                             signature = c("CD3E", "CD3D", 
                                           "CD2", "GZMK"), 
                             negative = T,
                             positive = F,
                             level = 1
)
scGate_model <- gating_model(model = scGate_model,
                             name = "Plasma", 
                             signature = c("DERL3", "ITM2C"), 
                             negative = T,
                             positive = F,
                             level = 2
)

scGate_model <- gating_model(model = scGate_model,
                             name = "Bcells", 
                             signature = c("CD79A", "CD79B", 
                                           "MS4A1", "BANK1"),
                             positive = T,
                             negative = F,                             
                             level = 3
)


merged.roider.singlets <- scGate(data = merged.roider.singlets, 
                        model = scGate_model, 
                        pos.thr = 0.25,
                        neg.thr = 0.15,
                        verbose = T,
                        ncores = 4,
                        pca.dim = 40,
                        reduction = "umap.unintegrated")
Computing UCell scores for all signatures using RNA assay...
# Visualize "Pure" and "Impure" cells
DimPlot(merged.roider.singlets, group.by = "is.pure")

# Subset and remove non-B cells
merged.roider.singlets <-  subset(merged.roider.singlets, is.pure %in% c("Pure"))

dim(merged.roider.singlets)
[1] 45068 27981

Remove cells with very high nCount_RNA values, set other final QC filters

merged.roider.singlets <- subset(merged.roider.singlets,
                                subset = nCount_RNA < 40000 & nCount_RNA > 500 & nFeature_RNA > 250)
dim(merged.roider.singlets)
[1] 45068 27894

Re-compute PCA

Re-scale data now that it has been subset

merged.roider.singlets[["RNA"]] <- split(merged.roider.singlets[["RNA"]], f = merged.roider.singlets$Sample)

merged.roider.singlets <- FindVariableFeatures(merged.roider.singlets, 
                                    assay="RNA", 
                                    layer="data", 
                                    selection.method = "vst", 
                                    nfeatures = 5000)
merged.roider.singlets <- ScaleData(merged.roider.singlets, vars.to.regress = "percent.mt")
merged.roider.singlets <- RunPCA(merged.roider.singlets, npcs = 200, verbose = FALSE)

Run Harmony integration

integrated.roider <- IntegrateLayers(merged.roider.singlets, 
                                method = HarmonyIntegration, 
                                orig.reduction = "pca",
                                new.reduction = "integrated.harmony"
)
Warning: HarmonyMatrix is deprecated and will be removed in the future from the
API in the future
Warning: Warning: The parameters do_pca and npcs are deprecated. They will be ignored for this function call and please remove parameters do_pca and npcs and pass to harmony cell_embeddings directly.
This warning is displayed once per session.
Warning: Warning: The parameter tau is deprecated. It will be ignored for this function call and please remove parameter tau in future function calls. Advanced users can set value of parameter tau by using parameter .options and function harmony_options().
This warning is displayed once per session.
Warning: Warning: The parameter block.size is deprecated. It will be ignored for this function call and please remove parameter block.size in future function calls. Advanced users can set value of parameter block.size by using parameter .options and function harmony_options().
This warning is displayed once per session.
Warning: Warning: The parameter max.iter.harmony is replaced with parameter max_iter. It will be ignored for this function call and please use parameter max_iter in future function calls.
This warning is displayed once per session.
Warning: Warning: The parameter max.iter.cluster is deprecated. It will be ignored for this function call and please remove parameter max.iter.cluster in future function calls. Advanced users can set value of parameter max.iter.cluster by using parameter .options and function harmony_options().
This warning is displayed once per session.
Warning: Warning: The parameter epsilon.cluster is deprecated. It will be ignored for this function call and please remove parameter epsilon.cluster in future function calls. Advanced users can set value of parameter epsilon.cluster by using parameter .options and function harmony_options().
This warning is displayed once per session.
Warning: Warning: The parameter epsilon.harmony is deprecated. It will be ignored for this function call and please remove parameter epsilon.harmony in future function calls. If users want to control if harmony would stop early or not, use parameter early_stop. Advanced users can set value of parameter epsilon.harmony by using parameter .options and function harmony_options().
This warning is displayed once per session.
Transposing data matrix
Using automatic lambda estimation
Initializing state using k-means centroids initialization
Harmony 1/10
Harmony 2/10
Harmony converged after 2 iterations

Save integrated object

saveRDS(integrated.roider,"Roider_DLBCL_scRNA_singlets_integrated_harmony.rds")

Identify clusters

integrated.roider <- FindNeighbors(integrated.roider, dims = 1:30, reduction = "integrated.harmony")
Computing nearest neighbor graph
Computing SNN
integrated.roider <- FindClusters(integrated.roider, 
                             resolution = seq(0.05,1,0.05), 
                             algorithm = 2)
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 27948
Number of edges: 873989

Running Louvain algorithm with multilevel refinement...
Maximum modularity in 10 random starts: 0.9688
Number of communities: 3
Elapsed time: 6 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 27948
Number of edges: 873989

Running Louvain algorithm with multilevel refinement...
Maximum modularity in 10 random starts: 0.9504
Number of communities: 4
Elapsed time: 6 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 27948
Number of edges: 873989

Running Louvain algorithm with multilevel refinement...
Maximum modularity in 10 random starts: 0.9368
Number of communities: 6
Elapsed time: 7 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 27948
Number of edges: 873989

Running Louvain algorithm with multilevel refinement...
Maximum modularity in 10 random starts: 0.9278
Number of communities: 7
Elapsed time: 6 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 27948
Number of edges: 873989

Running Louvain algorithm with multilevel refinement...
Maximum modularity in 10 random starts: 0.9193
Number of communities: 9
Elapsed time: 7 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 27948
Number of edges: 873989

Running Louvain algorithm with multilevel refinement...
Maximum modularity in 10 random starts: 0.9113
Number of communities: 9
Elapsed time: 7 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 27948
Number of edges: 873989

Running Louvain algorithm with multilevel refinement...
Maximum modularity in 10 random starts: 0.9034
Number of communities: 9
Elapsed time: 6 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 27948
Number of edges: 873989

Running Louvain algorithm with multilevel refinement...
Maximum modularity in 10 random starts: 0.8963
Number of communities: 11
Elapsed time: 6 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 27948
Number of edges: 873989

Running Louvain algorithm with multilevel refinement...
Maximum modularity in 10 random starts: 0.8895
Number of communities: 12
Elapsed time: 6 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 27948
Number of edges: 873989

Running Louvain algorithm with multilevel refinement...
Maximum modularity in 10 random starts: 0.8830
Number of communities: 12
Elapsed time: 8 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 27948
Number of edges: 873989

Running Louvain algorithm with multilevel refinement...
Maximum modularity in 10 random starts: 0.8766
Number of communities: 12
Elapsed time: 6 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 27948
Number of edges: 873989

Running Louvain algorithm with multilevel refinement...
Maximum modularity in 10 random starts: 0.8701
Number of communities: 12
Elapsed time: 7 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 27948
Number of edges: 873989

Running Louvain algorithm with multilevel refinement...
Maximum modularity in 10 random starts: 0.8638
Number of communities: 13
Elapsed time: 6 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 27948
Number of edges: 873989

Running Louvain algorithm with multilevel refinement...
Maximum modularity in 10 random starts: 0.8583
Number of communities: 13
Elapsed time: 6 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 27948
Number of edges: 873989

Running Louvain algorithm with multilevel refinement...
Maximum modularity in 10 random starts: 0.8531
Number of communities: 14
Elapsed time: 7 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 27948
Number of edges: 873989

Running Louvain algorithm with multilevel refinement...
Maximum modularity in 10 random starts: 0.8483
Number of communities: 15
Elapsed time: 7 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 27948
Number of edges: 873989

Running Louvain algorithm with multilevel refinement...
Maximum modularity in 10 random starts: 0.8432
Number of communities: 15
Elapsed time: 6 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 27948
Number of edges: 873989

Running Louvain algorithm with multilevel refinement...
Maximum modularity in 10 random starts: 0.8384
Number of communities: 16
Elapsed time: 6 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 27948
Number of edges: 873989

Running Louvain algorithm with multilevel refinement...
Maximum modularity in 10 random starts: 0.8336
Number of communities: 15
Elapsed time: 6 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 27948
Number of edges: 873989

Running Louvain algorithm with multilevel refinement...
Maximum modularity in 10 random starts: 0.8292
Number of communities: 17
Elapsed time: 7 seconds
integrated.roider <- RunUMAP(integrated.roider,
                        reduction = "integrated.harmony",
                        dims = 1:30,
                        reduction.name = "umap.harmony")
11:08:18 UMAP embedding parameters a = 0.9922 b = 1.112
Found more than one class "dist" in cache; using the first, from namespace 'BiocGenerics'
Also defined by 'spam'
11:08:18 Read 27948 rows and found 30 numeric columns
11:08:18 Using Annoy for neighbor search, n_neighbors = 30
Found more than one class "dist" in cache; using the first, from namespace 'BiocGenerics'
Also defined by 'spam'
11:08:18 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
11:08:26 Writing NN index file to temp file /tmp/RtmpGN9GPk/file150e28783fb059
11:08:26 Searching Annoy index using 1 thread, search_k = 3000
11:08:37 Annoy recall = 100%
11:08:39 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
11:08:43 Initializing from normalized Laplacian + noise (using RSpectra)
11:08:43 Commencing optimization for 200 epochs, with 1229528 positive edges
11:09:06 Optimization finished

Plot cluster results across range of resolutions

DimPlot(integrated.roider,
        reduction = "umap.harmony", 
        group.by = paste0("RNA_snn_res.",seq(0.05,1,0.05)), 
        label = TRUE)

Plot clusters and proceed using res = 0.5 as an example

Idents(integrated.roider) <- integrated.roider$RNA_snn_res.0.5
integrated.roider$seurat_clusters <- integrated.roider$RNA_snn_res.0.5
table(integrated.roider$seurat_clusters)

   0    1    2    3    4    5    6    7    8    9   10   11 
7348 6288 3436 2410 1906 1900 1869 1085  852  489  278   87 
DimPlot(integrated.roider,reduction = "umap.harmony", group.by = "RNA_snn_res.0.5", label = TRUE)

DimPlot(integrated.roider,reduction = "umap.harmony", group.by = "Disease")

DimPlot(integrated.roider,reduction = "umap.harmony", group.by = "Sample")

Plot QC

VlnPlot(integrated.roider, features = "nCount_RNA", group.by="seurat_clusters",pt.size=0) + NoLegend()

VlnPlot(integrated.roider, features = "nFeature_RNA", group.by="seurat_clusters",pt.size=0) + NoLegend()

VlnPlot(integrated.roider, features = "percent.mt", group.by="seurat_clusters",pt.size=0) + NoLegend()

Join with published cluster labels

Note their annotations are provided only on a subset of cells, about 6,000 out of 26,000 total B cells. This means most of our cells have NA for cluster here. We could request cluster labels for all cells from the authors to improve this

overlap <- intersect(rownames(roider.meta$bcells),colnames(integrated.roider))
length(overlap)
[1] 5565
roider.meta.subset <- roider.meta$bcells[overlap,]
pub.cluster <- roider.meta.subset[colnames(integrated.roider),3]
integrated.roider <- AddMetaData(integrated.roider, pub.cluster, col.name="roider_clusters")

DimPlot(integrated.roider,
    group.by = "roider_clusters",
    reduction = "umap.harmony",
    cells = overlap
    ) + 
DimPlot(integrated.roider,
    group.by = "seurat_clusters",
    reduction = "umap.harmony",
    cells = overlap
    )

Plot Roider’s UMAP coordinates for these overlapping cells but colored by our clusters

rownames_to_column(roider.meta.subset,var="Cell") %>%
    as_tibble() %>%
    inner_join(enframe(integrated.roider$seurat_clusters[overlap]) %>%
                   dplyr::rename("Cell" = name, "our.cluster" = value),
               by = "Cell"
               ) %>%
    ggplot(aes(x=UMAP1,y=UMAP2,color=our.cluster)) +
    geom_point(size=0.5) +
    theme_classic()

Make Sankey plot comparing cluster labels

roider.clusts <- enframe(integrated.roider$roider_clusters[overlap]) %>% 
    dplyr::rename("Barcode" = name, "Roider.cluster" = value) %>%
    mutate(Roider.cluster = as.numeric(as.character(Roider.cluster)))

our.clusts <- enframe(integrated.roider$seurat_clusters[overlap]) %>% 
    dplyr::rename("Barcode" = name, "our.cluster" = value) %>%
    mutate(our.cluster = as.numeric(as.character(our.cluster)))


make_long (
    full_join(roider.clusts,our.clusts,by="Barcode"),
    Roider.cluster,our.cluster) %>%
    ggplot(aes(x = x, 
        next_x = next_x, 
        node = node, 
        next_node = next_node,
        fill = factor(node),
        label = node)) +
        geom_sankey() +
        theme_sankey(base_size = 18) +
        theme(legend.position = "none") +
        xlab("") +
        geom_sankey_text(size = 3, color = "black") +
        scale_x_discrete(breaks = c("Roider.cluster","our.cluster"),labels = c("Roider clusters","Our clusters"))

Identify cursory marker genes of each cluster at this resolution

integrated.roider[['RNA']] <- JoinLayers(integrated.roider[['RNA']])
DefaultAssay(integrated.roider) <- "RNA"

vargenes <- presto::wilcoxauc(integrated.roider, 'seurat_clusters', seurat_assay = 'RNA')
top_vargenes <- top_markers(vargenes, n = 100, auc_min = 0.5, pct_in_min = 20, pct_out_max = 20)
top_vargenes
# A tibble: 100 × 13
    rank `0`   `1`   `10`  `11`  `2`   `3`   `4`   `5`   `6`   `7`   `8`   `9`  
   <int> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
 1     1 HSPA… IGLC7 MIR1… GNG11 ATP5… UBE2C CD9   SRGN  AC11… GPR1… FBLN5 AL92…
 2     2 HSPA… IGHE… EBI3  PCDH9 IGHG3 MKI67 RP5-… LAG3  SELE… TNFR… RP5-… SERP…
 3     3 CH17… AL92… IL4I1 GADD… ATP5… NUSA… RP11… SMIM… C1or… LINC… DNTT  MAP3…
 4     4 RP11… DMD   PLEK  FBP1  NOP53 CKS1B CCDC… S100… PLAC8 AC24… CCDC… GOLT…
 5     5 RP11… MIR3… CD44  CCDC… ATP5… ZWINT YBX3  ATP5… ATP5… CD44  CXor… TCL1B
 6     6 TNF   ITPR2 NFKB2 ZEB2  SELE… SMC4  KIAA… JSRP1 PLPP5 RIPO… ACY3  CD1C 
 7     7 RP11… EIF1… MARC… CSNK… SELE… BIRC5 PLAC8 IFI27 AC24… ATP5… CD9   IFNG…
 8     8 IGHGP LINC… SLAM… EML6  ELOB  TOP2A DNTT  MDK   AC00… ATP5… KIAA… TYMS 
 9     9 HIST… BFSP2 BATF  BMP3  ATP5… CENPF SOX4  CTSD  SLC2… AC11… PLAC8 HTR3A
10    10 SGK1  CCDC… ICAM1 RP11… ATP5… TYMS  RASS… LINC… RIPO… PLAC8 XIST  SCOC 
# ℹ 90 more rows
write_tsv(top_vargenes,"Roider_DLBCL_scRNA_singlets_integrated_harmony_clustered_prestoMarkers.tsv")

Plot top markers identified as a dotplot

top_vargenes <- top_markers(vargenes, n = 5, auc_min = 0.5, pct_in_min = 25, pct_out_max = 10)
top_markers <- top_vargenes %>%
        select(-rank) %>% 
        unclass() %>% 
        stack() %>%
        pull(values) %>%
        unique() %>%
        .[!is.na(.)]

# Add genes used in B-cell gating above
top_markers <- unique(c(top_markers,"CD79A","CD79B","MS4A1","BANK1","CD3E","CD2","CD3D","GZMK","DERL3", "ITM2C","CTLA4","TIGIT","LAG3","PRF1","SRGN"))

# Compute aggregated expression values of these genes and cluster them to order the figure
rna <- AverageExpression(integrated.roider,assay="RNA",slot="data")
As of Seurat v5, we recommend using AggregateExpression to perform pseudo-bulk analysis.
First group.by variable `ident` starts with a number, appending `g` to ensure valid variable names
This message is displayed once per session.
rna.sub <- rna$RNA[top_markers,]
cors.genes <- as.dist(1-cor(as.matrix(t(rna.sub)),method="pearson"))
hc.genes <- hclust(cors.genes)
top_markers.sorted <- rownames(rna.sub)[hc.genes$order]

# Plot
DotPlot(integrated.roider,features=top_markers.sorted,assay="RNA",cols=c("blue","red"),cluster.idents=T) + RotatedAxis()

Plot proportions per cluster per disease type

as.data.frame(table(integrated.roider$seurat_clusters, integrated.roider$Sample)) %>% 
  as_tibble() %>%
  dplyr::rename("Sample" = Var2,"Cluster" = Var1) %>%
  tidyr::extract(Sample,into=c("Disease","Patient"),regex="(.+)(\\d)",remove=F) %>%
  group_by(Sample) %>%
  mutate(Proportion = Freq / sum(Freq)) %>%
  ggplot(aes(fill = Disease, x = Disease, y=Proportion)) +
  geom_boxplot(outlier.shape=NA) +
  geom_point(pch = 21, position = position_jitterdodge(jitter.width=0.2)) +
  facet_wrap(~Cluster,scales="free")

Compute GCB module score per cell and plot in UMAP space and as dotplot

Derived from HCATonsilData v2 for GCB clusters

combinedGCB <- unique(c("AICDA","MME","RGS13","MEF2B","STMN1","MKI67","PCLAF","NUSAP1","TOP2A","HMGB2","SUGCT","MEF2B","NEIL1","SEC14L1","AICDA","BCL7A","MME","BACH2","TCL1A","CD38","MEF2B","RGS13","SERPINA9","LMO2","HMCES","KLHL6","LRMP","FGD6","MARCKSL1","NEIL1","LMO2","RGS13","MEF2B","BCL2A1","SERPINA9","GMDS","FGD6","RFTN1","PLEK","LRMP","BATF3","MYCBP","NME1","FABP5","PCLAF","MCM7","EBI3","PAICS","PCNA","LMO2"))
integrated.roider <- AddModuleScore(integrated.roider, features = list(combinedGCB), name = "GCB")

UMAP

FeaturePlot(integrated.roider,  reduction = "umap.harmony",  features = "GCB1", order=T) +
    FeaturePlot(integrated.roider,  reduction = "umap.harmony",  features = "HLA-E", order=T) +
    plot_layout(ncol=2)

Dotplot

DotPlot(integrated.roider,features=combinedGCB,assay="RNA",cols=c("blue","red"),cluster.idents=T) + RotatedAxis()

Average GCB scores by cluster and sort

as_tibble(rownames_to_column(integrated.roider@meta.data,var="Cells")) %>%
    dplyr::select(Cells,seurat_clusters,GCB1) %>%
    group_by(seurat_clusters) %>%
    summarize(GCB = mean(GCB1)) %>%
    dplyr::arrange(desc(GCB))
# A tibble: 12 × 2
   seurat_clusters      GCB
   <fct>              <dbl>
 1 3                0.289  
 2 4                0.117  
 3 8                0.103  
 4 9                0.0605 
 5 11               0.0451 
 6 2                0.0201 
 7 1                0.0102 
 8 10              -0.00964
 9 0               -0.0555 
10 6               -0.120  
11 5               -0.167  
12 7               -0.193  

Plot GCB signature vs HLA-E expression

Define which clusters were represented >5% in DLBCL tumors

DLBCLclusters <- as.data.frame(table(integrated.roider$seurat_clusters, integrated.roider$Sample)) %>% 
    as_tibble() %>%
    dplyr::rename("Sample" = Var2,"Cluster" = Var1) %>%
    tidyr::extract(Sample,into=c("Disease","Patient"),regex="(.+)(\\d)",remove=F) %>%
    group_by(Sample) %>%
    mutate(Proportion = Freq / sum(Freq)) %>% 
    group_by(Disease,Cluster) %>% 
    summarize(Mean = mean(Proportion)) %>% 
    dplyr::filter(Disease=="DLBCL") %>% 
    dplyr::filter(Mean > 0.05) %>% 
    pull(Cluster) %>% 
    as.character() %>% 
    as.numeric()
`summarise()` has grouped output by 'Disease'. You can override using the
`.groups` argument.
DLBCLclusters
[1] 0 1 3 4 5 8

Plot boxplot of HLA-E expression divided by GCB vs non-GCB cell clusters, each point is a sample

# Get top scoring cluster for GCB signature
gcb_clust <- as_tibble(rownames_to_column(integrated.roider@meta.data,var="Cells")) %>%
    dplyr::select(Cells,seurat_clusters,GCB1) %>%
    group_by(seurat_clusters) %>%
    summarize(GCB = mean(GCB1)) %>%
    slice_max(GCB, n = 1) %>%
    pull(seurat_clusters)
  
as_tibble(rownames_to_column(integrated.roider@meta.data,var="Cells")) %>%
    dplyr::select(Cells,Disease,Sample,seurat_clusters,GCB1) %>%
    group_by(seurat_clusters,Disease,Sample) %>%
    left_join(enframe(integrated.roider@assays$RNA$data["HLA-E",],name = "Cells",value = "HLAE") %>% as_tibble(), by="Cells") %>%
    summarize(GCB = mean(GCB1),HLAE = mean(HLAE)) %>%
    dplyr::filter(Disease=="DLBCL" & seurat_clusters %in% DLBCLclusters) %>%
    mutate(CellClass = case_when(
        seurat_clusters == gcb_clust ~ "GCB",
        T ~ "NonGCB"
        )
    ) %>%
    ggplot(aes(x=CellClass,y=HLAE,fill=CellClass)) + 
    geom_boxplot(outlier.shape=NA) +
    geom_point(pch = 21, position = position_jitterdodge(jitter.width=0.2))
`summarise()` has grouped output by 'seurat_clusters', 'Disease'. You can
override using the `.groups` argument.

Repeat above but plot GCB from rLN vs all other cells

as_tibble(rownames_to_column(integrated.roider@meta.data,var="Cells")) %>%
    dplyr::select(Cells,Disease,Sample,seurat_clusters,GCB1) %>%
    left_join(enframe(integrated.roider@assays$RNA$data["HLA-E",],name = "Cells",value = "HLAE") %>% as_tibble(), by="Cells") %>%
    mutate(CellClass = case_when(
        seurat_clusters == gcb_clust & Disease == "rLN" ~ "GCB",
        seurat_clusters == gcb_clust & Disease != "rLN" ~ "GCBother",
        T ~ "NonGCB"
      )
    ) %>%
    group_by(CellClass,Disease,Sample) %>%
    summarize(GCB = mean(GCB1),HLAE = mean(HLAE)) %>%
    ggplot(aes(x=Disease,y=HLAE,fill=Disease)) +
    geom_boxplot(outlier.shape=NA) +
    geom_point(pch = 21, position = position_jitterdodge(jitter.width=0.2)) +
    geom_text(aes(label=Sample),hjust=0.5,vjust=0.5) +
    facet_grid(~CellClass,scales="free_x",space = "free_x")
`summarise()` has grouped output by 'CellClass', 'Disease'. You can override
using the `.groups` argument.

Plot scatterplot of GCB signature vs HLA-E expression, colored by GCB vs non-GCB, each point is a cell

as_tibble(rownames_to_column(integrated.roider@meta.data,var="Cells")) %>%
    dplyr::select(Cells,Disease,Sample,seurat_clusters,GCB1) %>%
    group_by(seurat_clusters,Disease,Sample) %>%
    left_join(enframe(integrated.roider@assays$RNA$data["HLA-E",],name = "Cells",value = "HLAE") %>% as_tibble(), by="Cells") %>%
    dplyr::filter(Disease=="DLBCL" & seurat_clusters %in% DLBCLclusters) %>%
    mutate(CellClass = case_when(
        seurat_clusters == gcb_clust ~ "GCB",
        T ~ "NonGCB"
    )
    ) %>%
    ggplot(aes(x = GCB1, y = HLAE, color = CellClass)) + 
    geom_point(cex=0.5) +
    geom_smooth(method = "lm") +
    theme_classic() +
    xlab("GCB signature score") +
    ylab("HLA-E expression") +
    scale_color_manual(values=c("#B07AA1","#499894")) +
    stat_cor(aes(color = CellClass), method = "pearson", label.x = 0.2)
`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'
png 
  2 

Plot ridgeline showing HLA-E expression divided by GCB vs non-GCB

as_tibble(rownames_to_column(integrated.roider@meta.data,var="Cells")) %>%
    dplyr::select(Cells,Disease,Sample,seurat_clusters,GCB1) %>%
    group_by(seurat_clusters,Disease,Sample) %>%
    left_join(enframe(integrated.roider@assays$RNA$data["HLA-E",],name = "Cells",value = "HLAE") %>% as_tibble(), by="Cells") %>%
    dplyr::filter(Disease=="DLBCL" & seurat_clusters %in% DLBCLclusters) %>%
    mutate(CellClass = case_when(
        seurat_clusters == gcb_clust ~ "GCB",
        T ~ "NonGCB"
    )
    ) %>%
    ggplot(aes(x = HLAE, y = CellClass, fill = CellClass)) + 
    geom_density_ridges(alpha=0.7) + 
    scale_x_continuous(expand = c(0, 0)) +
    scale_y_discrete(expand = c(0, 0)) +
    coord_cartesian(clip = "off") + 
    theme_ridges(grid = FALSE, center_axis_labels = TRUE)
Picking joint bandwidth of 0.127

Save clustered object

saveRDS(integrated.roider,"Roider_DLBCL_scRNA_singlets_integrated_harmony_clustered.rds")

Get session info

sessionInfo()
R version 4.3.1 (2023-06-16)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Rocky Linux 8.10 (Green Obsidian)

Matrix products: default
BLAS/LAPACK: /usr/lib64/libopenblasp-r0.3.15.so;  LAPACK version 3.9.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

time zone: America/New_York
tzcode source: system (glibc)

attached base packages:
[1] stats4    stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] ggpubr_0.6.0                ggsankey_0.0.99999         
 [3] scGate_1.6.0                ggridges_0.5.6             
 [5] scDblFinder_1.14.0          harmony_1.2.0              
 [7] vctrs_0.6.5                 patchwork_1.3.0            
 [9] scater_1.34.0               scuttle_1.10.3             
[11] speckle_1.0.0               Matrix_1.6-4               
[13] readxl_1.4.3                SingleCellExperiment_1.22.0
[15] SummarizedExperiment_1.30.2 Biobase_2.60.0             
[17] GenomicRanges_1.52.1        GenomeInfoDb_1.36.4        
[19] IRanges_2.34.1              S4Vectors_0.38.2           
[21] BiocGenerics_0.46.0         MatrixGenerics_1.12.3      
[23] matrixStats_1.2.0           flexmix_2.3-19             
[25] lattice_0.22-5              SeuratWrappers_0.3.19      
[27] miQC_1.8.0                  lubridate_1.9.3            
[29] forcats_1.0.0               stringr_1.5.1              
[31] dplyr_1.1.4                 purrr_1.0.2                
[33] readr_2.1.5                 tidyr_1.3.1                
[35] tibble_3.2.1                ggplot2_3.5.1              
[37] tidyverse_2.0.0             Seurat_5.1.0               
[39] SeuratObject_5.0.2          sp_2.1-3                   
[41] sctransform_0.4.1           glmGamPoi_1.12.2           
[43] presto_1.0.0                Rcpp_1.0.12                
[45] devtools_2.4.5              usethis_2.2.2              
[47] data.table_1.15.0          

loaded via a namespace (and not attached):
  [1] fs_1.6.3                  spatstat.sparse_3.0-3    
  [3] bitops_1.0-7              httr_1.4.7               
  [5] RColorBrewer_1.1-3        backports_1.4.1          
  [7] profvis_0.3.8             tools_4.3.1              
  [9] utf8_1.2.4                R6_2.5.1                 
 [11] mgcv_1.9-1                lazyeval_0.2.2           
 [13] uwot_0.1.16               urlchecker_1.0.1         
 [15] withr_3.0.0               gridExtra_2.3            
 [17] progressr_0.14.0          cli_3.6.2                
 [19] spatstat.explore_3.2-6    fastDummies_1.7.3        
 [21] labeling_0.4.3            spatstat.data_3.0-4      
 [23] pbapply_1.7-2             Rsamtools_2.16.0         
 [25] R.utils_2.12.3            parallelly_1.37.0        
 [27] sessioninfo_1.2.2         limma_3.56.2             
 [29] generics_0.1.3            BiocIO_1.10.0            
 [31] vroom_1.6.5               ica_1.0-3                
 [33] spatstat.random_3.2-2     car_3.1-2                
 [35] ggbeeswarm_0.7.2          fansi_1.0.6              
 [37] abind_1.4-5               R.methodsS3_1.8.2        
 [39] lifecycle_1.0.4           yaml_2.3.8               
 [41] edgeR_3.42.4              carData_3.0-5            
 [43] SparseArray_1.2.2         Rtsne_0.17               
 [45] grid_4.3.1                promises_1.2.1           
 [47] dqrng_0.3.2               crayon_1.5.2             
 [49] miniUI_0.1.1.1            beachmat_2.16.0          
 [51] cowplot_1.1.3             metapod_1.8.0            
 [53] pillar_1.9.0              knitr_1.45               
 [55] rjson_0.2.21              xgboost_1.7.7.1          
 [57] future.apply_1.11.1       codetools_0.2-19         
 [59] leiden_0.4.3.1            glue_1.7.0               
 [61] remotes_2.4.2.1           png_0.1-8                
 [63] spam_2.10-0               cellranger_1.1.0         
 [65] gtable_0.3.4              cachem_1.0.8             
 [67] xfun_0.42                 S4Arrays_1.2.0           
 [69] mime_0.12                 survival_3.5-8           
 [71] statmod_1.5.0             bluster_1.10.0           
 [73] ellipsis_0.3.2            fitdistrplus_1.1-11      
 [75] ROCR_1.0-11               nlme_3.1-164             
 [77] bit64_4.0.5               RcppAnnoy_0.0.22         
 [79] irlba_2.3.5.1             vipor_0.4.7              
 [81] KernSmooth_2.23-22        colorspace_2.1-0         
 [83] nnet_7.3-19               ggrastr_1.0.2            
 [85] UCell_2.6.2               tidyselect_1.2.0         
 [87] bit_4.0.5                 compiler_4.3.1           
 [89] BiocNeighbors_1.18.0      DelayedArray_0.26.7      
 [91] plotly_4.10.4             rtracklayer_1.60.1       
 [93] scales_1.3.0              lmtest_0.9-40            
 [95] digest_0.6.34             goftest_1.2-3            
 [97] spatstat.utils_3.0-4      rmarkdown_2.25           
 [99] XVector_0.40.0            htmltools_0.5.7          
[101] pkgconfig_2.0.3           sparseMatrixStats_1.12.2 
[103] fastmap_1.1.1             rlang_1.1.3              
[105] htmlwidgets_1.6.4         shiny_1.8.0              
[107] DelayedMatrixStats_1.22.6 farver_2.1.1             
[109] zoo_1.8-12                jsonlite_1.8.8           
[111] BiocParallel_1.34.2       R.oo_1.26.0              
[113] BiocSingular_1.16.0       RCurl_1.98-1.14          
[115] magrittr_2.0.3            modeltools_0.2-23        
[117] GenomeInfoDbData_1.2.10   dotCall64_1.1-1          
[119] munsell_0.5.0             viridis_0.6.5            
[121] reticulate_1.35.0         stringi_1.8.3            
[123] zlibbioc_1.46.0           MASS_7.3-60.0.1          
[125] plyr_1.8.9                pkgbuild_1.4.3           
[127] parallel_4.3.1            listenv_0.9.1            
[129] ggrepel_0.9.5             deldir_2.0-2             
[131] Biostrings_2.68.1         splines_4.3.1            
[133] tensor_1.5                hms_1.1.3                
[135] locfit_1.5-9.8            igraph_2.0.2             
[137] spatstat.geom_3.2-8       ggsignif_0.6.4           
[139] RcppHNSW_0.6.0            reshape2_1.4.4           
[141] ScaledMatrix_1.8.1        pkgload_1.3.4            
[143] XML_3.99-0.16.1           evaluate_0.23            
[145] scran_1.28.2              BiocManager_1.30.22      
[147] tzdb_0.4.0                httpuv_1.6.14            
[149] RANN_2.6.1                polyclip_1.10-6          
[151] future_1.33.1             scattermore_1.2          
[153] rsvd_1.0.5                broom_1.0.5              
[155] xtable_1.8-4              restfulr_0.0.15          
[157] RSpectra_0.16-1           rstatix_0.7.2            
[159] later_1.3.2               viridisLite_0.4.2        
[161] memoise_2.0.1             beeswarm_0.4.0           
[163] GenomicAlignments_1.36.0  cluster_2.1.6            
[165] timechange_0.3.0          globals_0.16.2